#include <stdexcept>
#include "../Iteration_Params.hxx"

Iteration_Params::Iteration_Params(const double &Xi, const double &Et,
                                   const double &Q, const bool &kxi,
                                   const bool &ket,
                                   const double &sin_dip,
                                   const double &cos_dip):
  xi(Xi), et(Et), q(Q), xi2(xi*xi), et2(et*et), q2(q*q), r2(xi2+et2+q2),
  r(sqrt(r2)), r3(r*r2), r5(r3*r2), y(et*cos_dip + q*sin_dip),
  d(et*sin_dip - q*cos_dip), tt(q==0 ? 0 : atan(xi*et/(q*r)))
{
  if(r==0)
    throw std::runtime_error("Singular radius in Iteration_Params");

  if(kxi)
    {
      alx=-log(r-xi);
      x11=0;
      x32=0;
    }
  else
    {
      double rxi=r+xi;
      alx=log(rxi);
      x11=1/(r*rxi);
      x32=(r+rxi)*x11*x11/r;
    }

  if(ket)
    {
      ale=-log(r-et);
      y11=0;
      y32=0;
    }
  else
    {
      double ret=r+et;
      ale=log(ret);
      y11=1/(r*ret);
      y32=(r+ret)*y11*y11/r;
    }
  ey=sin_dip/r - y*q/r3;
  ez=cos_dip/r + d*q/r3;
  fy=d/r3 + xi2*y32*sin_dip;
  fz=y/r3 + xi2*y32*cos_dip;

  gy=2*x11*sin_dip - y*q*x32;
  gz=2*x11*cos_dip + d*q*x32;
  hy=d*q*x32 + xi*q*y32*sin_dip;
  hz=y*q*x32 + xi*q*y32*cos_dip;
}
